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Motivated by recent experimental progress in preparing encapsulated graphene sheets with ultra- 
high mobilities up to room temperature, we present a theoretical study of dc transport in doped 
graphene in the hydrodynamic regime. By using the continuity and Navier-Stokes equations, we 
demonstrate analytically that measurements of non-local resistances in multi-terminal Hall bar 
devices can be used to extract the hydrodynamic shear viscosity of the two-dimensional (2D) electron 
liquid in graphene. We also discuss how to probe the viscosity-dominated hydrodynamic transport 
regime by scanning probe potentiometry and magnetometry. Our approach enables measurements 
of the viscosity of any 2D electron liquid in the hydrodynamic transport regime. 


I. INTRODUCTION 

Transport in systems with many particles (such as 
gases and liquids) undergoing very frequent inter-particle 
collisions has been studied for more than two centuries 
and is described by the theory of hydrodynamics^^^. In 
the hydrodynamic regime, transport is described by^^^ 
three non-linear partial-differential equations—the con¬ 
tinuity, Navier-Stokes, and energy-transport equations— 
reflecting the conservation of mass, momentum and en¬ 
ergy, respectively. The Navier-Stokes equation contains 
two transport coeflicients^^^: the shear viscosity, rj, which 
describes friction between adjacent layers of fluid moving 
with different velocities, and the bulk viscosity, /, which 
describes dissipation arising in a liquid due to homoge¬ 
neous compression-like deformations. The energy trans¬ 
port equation contains the thermal conductivity k, which 
describes dissipative heat flow between regions with dif¬ 
ferent temperatures. These coefficients quantify the ten¬ 
dency of the liquid to restore a homogeneous state in 
response to a velocity or thermal gradient: they there¬ 
fore control the magnitude of non-local contributions to 
the linear-response functions of the liquid. 

Viscous flow of dilute classical gases attracted the at¬ 
tention of Maxwell, who theoretically discovered^ a puz¬ 
zling property of the shear viscosity of a dilute gas. Using 
a molecular approach^, he found that the shear viscosity 
of a dilute gas is independent of density n (and depends 
on temperature T according to 77 oc T^/^), a counterintu¬ 
itive result that he felt needed immediate experimental 
testing"^. The importance of rj in the hydrodynamic be¬ 
havior of dilute gases and liquids stems from the fact that 
this parameter controls the term of the Navier-Stokes 
equation that opposes turbulent flow^“^. 

Recent years have witnessed a tremendous interdisci¬ 
plinary interest in the hydrodynamic flow of strongly in¬ 
teracting quantum fluids. This interest was sparked by 
a series of results^, which were obtained via the anti- 
de Sitter/conformal held theory (AdS/CFT) correspon¬ 
dence, for the shear viscosity of a large class of strongly 


interacting thermal quantum held theories. These ef¬ 
forts culminated in 2005 when it was conjectured® that 
all quantum huids obey the following universal lower 
bound: ? 7 /s > ?i/(47rfcB)j where s is the entropy density. 
Note that this bound does not contain the speed of light, 
thereby explaining why the conjecture was extended also 
to non-relativistic quantum held theories. Fluids that 
saturate this bound have been dubbed “nearly perfect 
huids” (NPFs) ' , i.e. huids that dissipate the smallest pos¬ 
sible amount of energy and satisfy the laws of hydrody¬ 
namics at distances as short as the inter-particle spacing. 
Currently, two laboratory systems come closest to satu¬ 
rating the AdS/CFT bound: i) the quark-gluon plasma®, 
which is created at Brookhaven’s Relativistic Heavy Ion 
Collider and at CERN’s Large Hadron Collider by bash¬ 
ing heavy (e.g. gold and lead) ions together and ii) ul¬ 
tracold atomic Fermi gases®’^® (such as ®Li) close to a 
Feshbach resonance. Although mathematical counterex¬ 
amples have appeared in the literature^^, there are no 
known experimental violations of the AdS/CFT bound. 

The present work is motivated by the following ques¬ 
tions: Do electron liquids display hydrodynamic behav¬ 
ior? If so, how can it be experimentally proven that the 
electron system has entered the hydrodynamic transport 
regime? Once in the hydrodynamic regime, how can the 
shear viscosity of an electron liquid be measured in a 
solid-state device? Can an electron liquid in a solid-state 
device be a NPF? 

Hydrodynamics has been used for a long time to de¬ 
scribe transport of electrons in solid-state devices^^^^®. 
However, since the (Bloch) momentum of an electron 
in a solid is a poorly conserved quantity due to colli¬ 
sions against impurities, phonons, and structural defects 
in the crystal, experimental signatures of hydrodynamic 
electron flow are expected only in ultra-clean crystals, at 
sufficiently low temperatures. Second, electron-electron 
(e-e) interactions need to be sufficiently strong to ensure 
that the mean free path ^ee for e-e collisions is the short¬ 
est length scale in the problem, i.e. iee ^ £,W,vf/lo. 
Here, wp is the Fermi velocity, £ is the mean free path for 
momentum-non-conserving collisions, W is the sample 
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size, and uj the frequency of the external perturbation. 

Unfortunately, the low-temperature requirement 
(needed to mitigate the breakdown of momentum 
conservation in a solid-state environment) and the 
strong e-e interaction requirement are conflicting: at low 
temperatures, where £ is large, the mean free path £ee 
for e-e collisions is also very large due to Pauli blocking. 
Indeed, normal Fermi liquids at low temperatures^®^^® 
display very large values of £ee, i-e. £ee oc (Tp/T)^ for 
temperatures T Tp, with Tp the Fermi temperature. 
(In two spatial dimensions there is a very well known^® 
logarithmic correction that has been dropped.) These se¬ 
vere restrictions, imposed by a rigid Fermi surface on the 
phase space for two-body e-e collisions, can be relieved 
by increasing temperature. Indeed, ige quickly decreases 
for increasing T. Short e-e mean free paths therefore 
require to operate at sufficiently elevated temperatures. 
At such temperatures, strong scattering of electrons 
against optical phonons (e.g. in polar cystals such as 
GaAs) often leads to the unwanted inequality £ < £ee- 
Realizing hydrodynamic flow at “high” temperatures 
therefore requires not only ultra-clean crystals but also 
crystals where electron-phonon coupling is extremely 
weak. We note that it is in principle easier to reach the 
hydrodynamic transport regime in non-polar crystals 
such as graphene, where the dominating mechanism 
at large temperatures is scattering of electrons against 
acoustic phonons. In this case £ decays as I/T, i.e. slower 
than £ee. 

For these reasons, evidence of hydrodynamic trans¬ 
port in solid-state devices is, to the best of our knowl¬ 
edge, limited to early work carried out by Molenkamp 
and de Jong®®’®^ in electrostatically defined wires in 
GaAs/AlGaAs heterostructures. These authors mea¬ 
sured a non-monotonic dependence of the four-point lon¬ 
gitudinal resistivity pxx on the electronic temperature 
T, which is increased above the lattice temperature by 
using a large dc heating current. The decrease of pxx 
with increasing temperature above a certain value of 
T was attributed to e-e collisions®®’®^. This is the so- 
called “Gurzhi effect” and will be discussed extensively 
in this Article. More recently, indirect evidence of hydro- 
dynamic flow comes from an explanation®^ of Coulomb 
drag between two neutral graphene sheets®®, which dif¬ 
fers from that offered by the authors of Ref. 33. 

Recent experimental progress®^^®®, however, has made 
it possible to fabricate samples with ultra-high carrier 
transport mobilities up to room temperature. These are 
graphene sheets encapsulated between thin hexagonal 
boron nitride (hBN) slabs, which display ultra-high mo¬ 
bilities reaching 10® cm^/(Vs) in a wide range of tem¬ 
peratures up to room temperature. These values can be 
achieved in GaAs/AlGaAs heterostructures only below 
100 K because of polar phonon scattering^®. In addition, 
the finite electron mass and moderate doping required 
to achieve high mobilities limits the Fermi energy to val¬ 
ues below a few tens of meV, which makes the Fermi 
level smearing important even at liquid nitrogen temper¬ 


atures. Encapsulated samples®®^®® are special because 
electrons roaming in graphene suffer very weak scatter¬ 
ing against acoustic phonons^®^'*"^ and because hBN pro¬ 
vides an exceptionally clean and flat dielectric environ¬ 
ment for graphene^®. Furthermore, microscopic calcu¬ 
lations based on many-body diagrammatic perturbation 
theory'®®^®^® indicate that the e-e mean free path £ee in 
graphene is shorter than 400 nm in a wide range of car¬ 
rier densities and temperatures T > 150 K. We therefore 
conclude that hBN/graphene/hBN stacks are ideal sam¬ 
ples where the long-sought hydrodynamic regime can be 
unveiled and explored. Indeed, recent non-local trans¬ 
port measurements®’® carried out in high-quality encap¬ 
sulated single-layer (SLG) and bilayer (BLG) graphene 
samples have demonstrated that this is the case. The 
authors of Ref. 49 have reported evidence of hydrody¬ 
namic transport, showing that doped graphene exhibits 
an anomalous (negative) voltage drop near current injec¬ 
tion points, which has been attributed to the formation 
of whirlpools in the electron flow. From measurements 
of non-local signals, Bandurin et extracted the vis¬ 
cosity of graphene’s electron liquid and found it to be in 
quantitative agreement with many-body theory calcula¬ 
tions®®. 

In this Article, we present a fully-analytical theoretical 
study of non-local dc transport in the two-dimensional 
(2D) electron liquid in a graphene sheet in the hydro- 
dynamic regime. In Sect. II we present the theoretical 
framework that was used in Ref. 49 to interpret the ex¬ 
perimental results, i.e. a linearized steady-state hydro- 
dynamic approach based on the continuity and Navier- 
Stokes equations. Suitable boundary conditions for these 
hydrodynamic equations are discussed in Sect. II A. In 
Sect. Ill we present analytical solutions for longitudinal 
transport in a rectangular Hall bar and we discuss the 
dependence of the solutions on the boundary conditions, 
providing details on the Gurzhi effect. In Sect. IV we 
present analytical results for the spatial depedence of the 
2D electrical potential, non-local resistance, and current- 
induced magnetic field, as obtained by solving the hydro- 
dynamic equations with the free-surface boundary con¬ 
ditions, which we believe to be the appropriate bound¬ 
ary conditions for the linear-response regime. Finally, in 
Sect. V we present a summary of our main results and 
offer some perspectives. 

We remark that, in the present work, we focus only 
on doped SLG and BLG sheets, where the applicabil¬ 
ity of the Fermi-liquid theory is granted. However, it is 
believed that the hydrodynamic behavior of the semimet¬ 
als is particularly interesting®®’®®" when these are in the 
charge neutral state. In this case the Fermi surface 
shrinks to a point and Fermi liquid theory is not ap¬ 
plicable. For example, the authors of Ref. 17 have found 
that the ratio 77 /s for the 2D electron liquid in a graphene 
sheet at the charge neutrality point comes close to satu¬ 
rating the AdS/CFT bound. In nearly neutral semimet¬ 
als £ee is also short due to frequent collisions between 
thermally excited carriers (T 3> Tp). It is, however, very 
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well known®^ that any theory at the charge neutrality 
point must take into account the spatially inhomogeneous 
pattern of electron-hole puddles created by disorder^^. In 
the regime of sizeable doping we consider, we can safely 
ignore puddles in encapsulated graphene devices®^’^^. 

II. LINEARIZED STEADY-STATE 
HYDRODYNAMIC THEORY 

We consider a two-dimensional electron liquid in a 
doped SLG or BLG sheet, deep in the hydrodynamic 
transport regime (4e 'C i,W). For the sake of defi¬ 
niteness, we consider the Hall bar geometry sketched in 
Fig. 1. Since the energy-momentum dispersion of elec¬ 
trons in these systems is particle-hole symmetric®^, we 
assume, without loss of generality, that the sample hosts 
a back gate-controlled equilibrium electron density equal 
to h. (The charge density is —eh, —e being the elec¬ 
tron charge.) We neglect thermally-excited carriers and 
coupling between charge and heat flow®®, which is strong 
only at the charge neutrality point. Finally, we consider 
the linear response regime and steady-state transport. 

In this framework of approximations, the hydrody¬ 
namic transport equations^®’^^ for the 2D electron liquid 
greatly simplify and reduce to 

V-J(r) = 0, (1) 

and 

— V<^(r)-t ^VV(r) = . (2) 

m T 

In Eqs. (I) and (2) we have introduced the linearized 
steady-state particle current density J(r) = nv{r), 
where v{r) is the linearized steady-state fluid-element 
velocity. 

Eq. (1) is the continuity equation, while Eq. (2) is the 
Navier-Stokes equation. The latter contains three forces 
acting on a fluid element: i) the electric force —eE{r) = 
eV(/)(r), written in terms of the electric potential ^(r) 
in the 2D plane where electrons are moving, which is 
generated by the steady-state charge distribution n{r) 
in response to the drive current I; ii) the internal force 
due to the shear viscosity rj = r]{n, T) of the 2D electron 
liquid, here written in terms of the kinematic viscosity^^® 


and iii) friction exerted on a fluid element by agents ex¬ 
ternal to the electron liquid such as phonons and impu¬ 
rities, which dissipate the fluid-element momentum at a 
rate = l/T(h,T). The latter is a phenomenological 
parameter, which depends on h and T and is commonly 
used in modelling transport in semiconductor devices®®. 

In Eqs. (2) and (3) m is a suitable effective mass de¬ 
fined by: 



FIG. 1. (Color online) Schematic representation of the non¬ 
local transport setnp analyzed in this work. A dc cnrrent I 
is injected (red arrow) into an encapsnlated graphene Hall 
bar of width W. Current injection occnrs at a lateral con¬ 
tact located at x = xo and y = —TT/2. The same current 
is drained (blue arrow) at a contact located at x = —xq and 
y = —W12. Measurements of voltage drops AV “near” the 
current injection region are sensitive to the kinematic viscos¬ 
ity v of the two-dimensional massless Dirac fermion liquid. 
The notion of “vicinity” between voltage probe and current 
injector is defined by a crucial length scale, i.e. the vorticity 
diffusion length = ydd". Here r (exceeding 1 ps in high- 
quality encapsulated devices) represents a phenomenological 
scattering time due to momentum-non-conserving collisions 
of a fluid element (and not of single electrons). 


where rric = tik-p /up is the 2D massless Dirac fermion cy¬ 
clotron mass®^, fcp = being the Fermi wave number 
and Up ^ 10® m/s the Fermi velocity, and me is the bare 
electron mass in vacuum. 

Multiplying both members of Eq. (2) by r, we obtain 

— V(/)(r) + DlV^J(r) = J{r) . (5) 

e 

In Eq. (5) we have introduced the following characteristic 
length scale of the problem: 

Di, = y/vT . ( 6 ) 

For T = 1 ps (as in high-quality hBN/graphene/hBN 
samples) and zz = 0.1 m^/s (see Ref. 48) we obtain D^, « 
0.3 /rm. 

The physical significance of D^, can be understood as 
follows. We hrst note that we can rewrite V‘^J{r) by 
using the following identity: 

VV(r)=V[V-J(r)]-Vx [Vx J(r)] . (7) 

Because of (1), we can drop the first term on the right- 
hand side of Eq. (7). The second term is finite and related 
to the vorticity®’^ 

a;(r) = x J(r) = u){r)z 
n 


m = 


rric, for SLG, 

0.03 me, for BLG 


( 4 ) 


( 8 ) 



4 


which in 2D is oriented along the z axis. We can then 
rewrite Eq. (5) as following 

—V(j>(r) — nD^V x wfr) = J(r) . (9) 

e 

Taking the curl of Eq. (9) and using the identities V x 
V(/>(r) = 0 and V • a;(r) = 0 (the latter being valid 
because u;’s only non-vanishing component is along z, 
while V acts only on the 2D x-y plane), we obtain a 
damped-diffusion equation for the vorticity 

DlV^uj{r) = u}{r) . (10) 

We therefore see that Di, plays the role of a diffusion 
length for uj{r). 

In Eq. (5) we have also introduced a “Drude-like” con¬ 
ductivity, 


Since we are in the hydrodynamic regime, cto should not 
be confused^^ with the ordinary dc conductivity in the 
diffusive transport regime: once again, r = r(h, T) repre¬ 
sents a phenomenological parameter that should be fit to 
experimental data, as we discuss below in Sect. III. This 
naive description of momentum-non-conserving collisions 
in the hydrodynamic transport regime can be relaxed by 
following similar arguments to those in Ref. 19: this is 
however well beyond the scope of the present Article and 
will be the topic of future studies. In the absence of vis¬ 
cosity, Eq. (2) reduces to a local version of Ohm’s law, 
i.e. J(r) = (j){r)/e. 

Finally, we note that taking the divergence of Eq. (5) 
and making use of Eq. (1) we obtain the Laplace equa¬ 
tion V^(/)(r) = 0 for the electric potential 4'{r) on the 
2D plane. This should not be confused with the usual 
three-dimensional (3D) Poisson equation for the 3D elec¬ 
trostatic potential $(r,z), 

-I-$(r, z) = 47ren(r)(5(z) . (12) 

The 2D potential in Eq. (2) is (j>{r) = 4>(r,z = 0). On 
the right-hand side of Eq. (12) we note the steady-state 
charge density distribution —en(r) which occurs in the 
sample in response to the drive current I. Eq. (12) 
needs to be solved in 3D space with suitable bound¬ 
ary conditions—depending on the dielectric environment, 
gates, etc. surrounding the graphene sheet—if one is in¬ 
terested in determining n(r). In this Article we will focus 
our attention on J{r) and <p{r). 

Eqs. (1) and (5) will be used to describe transport in 
the Hall bar geometry pictorially represented in Fig. 1. 
Mathematically, it is convenient to work in a Hall bar 
of infinite length in the longitudinal direction x, since 
this allows us to use the Fourier transform to solve the 
equations of motion—see Sect. IV. The width W of the 
Hall bar will be kept finite. In the next Section we will 
describe a crucially important ingredient of the theory: 
boundary conditions. 


A. Boundary conditions 

In order to find (j){r) and J{r) in the Hall bar geome¬ 
try depicted in Fig. 1, we need to solve Eqs. (1) and (5) 
in the rectangle (— 00 , 00 ) x [—IV/2, IE/2], with appro¬ 
priate boundary conditions (BCs) at the edges, i.e. at 
y = ±IE/2. 

Lateral electrodes acting as current injectors/collectors 
are described through BCs on the component of the cur¬ 
rent perpendicular to the edges: 

Jyix,y = ±W/2) = J±{x). (13) 

Here J7± {x) is a function that describes a distribution of 
current injectors and collectors on the upper (lower) edge 
of the multi-terminal Hall bar. It is through Eq. (13) that 
the total drive current / injected into the system at the 
boundaries enters the problem. 

Following Abanin et , we model the electrodes as 
point-like (i.e. delta-function) sources and sinks. (A more 
realistic modeling of electrodes has been carried out in 
Ref. 49, where finite-width effects and metallic boundary 
conditions at extended electrodes have been taken into 
account in a fully numerical solution of Eqs. (1) and (5). 
Such details have essentially no impact on the physics we 
are going to highlight below.) For example, for the setup 
depicted in Fig. 1 with a current injector at a; = Xq, a 
current collector at x = —xq, and no injectors/collectors 
on the upper edge, we will use: 

J-{x) = --S{x - xo) +-S{x + xo) , (14) 

e e 

and J7+(x) = 0. 

In the presence of a finite shear viscosity ly, we need an 
additional BC on the tangential component of the current 
at the top {y = +Wl2) and bottom (y = —IE/2) edges 
of the Hall bar. We use the following BC: 

[dyUx,y)+d.Jyix,y)]^^^^^, = ^ Mx,y = ±W/2) 

(15) 

where 1^ is a “boundary slip length”, i.e. a length scale 
describing friction at the physical boundaries of the sam¬ 
ple. This BC can be explained as following. The left- 
hand side of Eq. (15) is proportional to the off-diagonal 
component of the stress tensor^, calculated at the edges 
of the Hall bar. It represents the tangential component 
of the frictional force exerted by the boundaries of the 
Hall bar on the 2D electron liquid^. This force depends 
on the tangential velocity of the 2D electron liquid and 
boundary roughness: in the linear-response regime, it is 
natural to replace such unknown dependence with a lin¬ 
ear law characterized by the single parameter 1^, as in 
the right-hand side of Eq. (15). 

In the description of transport of molecular liquids in 
constrained geometries, like water in a pipe, where the 
interactions between the molecules of the fluid and the 
walls of the container are of the same nature as of those 
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between molecules of the fluid, the most used BCs are 
the so-called “no-slip” BCs^, in which the component of 
the current tangential to the boundary vanishes. The 
no-slip BCs can be obtained from Eq. (15) by taking 
the limit Zb 0. In the opposite limit of a free-surface 
geometry, like the surface of water in an open bucket, the 
tangential force applied from the boundary to the fluid 
element vanishes at the boundary. These “free-surface” 
BCs^’^^ can be obtained from Eq. (15) by taking the limit 
Zb —^ -kOO. 

Which of these BCs should be used to model the exper¬ 
iments in Ref. 49 will become clear at the end of Sect. III. 


where we have introduced the length 

We can calculate the total longitudinal current / carried 
by the flow by integrating Eq. (17) in the transverse di¬ 
rection, i.e. 

I=-e dyUy) = aoWE,il-T) , (19) 

J-W/2 

where we have defined the dimensionless quantity 


B. Applicability of the linearized theory 


£>2 

.F.2^sinh 


W \ 

■ 


( 20 ) 


The validity of the linearized Navier-Stokes equation 
(5) relies on the smallness of the Reynolds number^^^ 
TZw- This is a dimensionless parameter (which depends 
on the sample geometry) that controls the smallness of 
the non-linear term [t;(r,t) • V]n(r,Z) in the convective 
derivative with respect to the viscous term. In our case 
we can define the Reynolds number as following: 


[v(r,t) ■ Vjv(r,t) 
iyV^v(r, t) 


vW 

ly 



env 


(16) 


where v is the typical value of the fluid-element velocity. 
For an injected current^® / = 2 x 10“^ A, a Hall bar width 
W = 1 ^m, and an equilibrium density h = 10^^ cm“^, 
we obtain v ^ I/[enW) « 10^ cm/s. We note that v 
is much smaller than the graphene Fermi velocity up 
10® m/s and the flow is therefore “non-relativistic.” The 
corresponding value of the Reynolds number is TZw 
10“® <C 1, obtained by using a kinematic viscosity v ~ 
10® cm^/s of the 2D electron liquid in graphene"^®. Our 
linearized theory in Eqs. (1) and (5) is therefore fully 
justified. 


III. LONGITUDINAL TRANSPORT AND THE 
GURZHI EFFECT 

We first consider the situation in which no current 
is injected or extracted laterally at the Hall bar edges, 
i.e. J±{x) = 0. 

In this case the local current J(r) does not depend on 
the longitudinal coordinate x and all the spatial deriva¬ 
tives with respect to x in Eqs. (1), (5), and (15) van¬ 
ish. The continuity equation implies that Jy does not 
depend on y and vanishes identically because of Eq. (13). 
Therefore also the y component of the electric field must 
vanish. The x component of the current respects the fol¬ 
lowing equation: Jx{y) - DldyJx{y) = -a^E^je, where 
E = — V(/(r) is the electric field. Note that E^ cannot 
depend on y because Ey vanishes and V x EZ = 0. The 
solution of this equation that fulfils the BC (15) is 


Measuring the longitudinal potential drop AV between 
two lateral contacts at positions x and x L yields a 
four-point longitudinal conductivity axx of the form: 

=CTo(l-J") . (21) 

Eq. (21) is the most important result of this Section. 

In the limit Zb —?► oo (i.e. free-surface BCs) —>■ 0. 
For this choice of BCs the longitudinal conductivity a^x 
depends only on the rate of momentum-non-conserving 
collisions (through gq) and is independent of v. 

On the other hand, in the limit Zb —> 0 (i.e. no-slip 
BCs) Eq. (21) reduces to 



^ W Y 

1 — 2 —- tanh 

[ W 

v2dJJ 


We can easily understand two asymptotic limits of 
Eq. (22). In the limit Dy <^W Eq. (22) yields Gxx = 
CTo(l — 2Dy/W): the small correction to the Drude-like 
conductivity gq is due to a reduction of the fluid-element 
velocity in a thin region of width Dy near the top and bot¬ 
tom edges of the Hall bar. In the opposite limit, Dy ^ 
W, we obtain Gxx = /{1201) — nW"^/{12mv). 

In this limit the problem is equivalent to that of Poiseuille 
flow in a pipe^’^, with a velocity profile Vx{y) that de¬ 
pends quadratically on the transverse coordinate y and 
a resistance that is entirely due to viscosity. 

A summary of our main results for longitudinal elec¬ 
tron transport in the presence of a Hnite viscosity is re¬ 
ported in Fig. 2. 


A. The Gurzhi effect 

We now would like to make a remark on the temper¬ 
ature dependence of Gxx in Eq. (21). For the sake of 
simplicity, we assume that Zb does not depend on tem¬ 
perature. We observe that the derivative of Gxx with 
respect to T, 


Jx{y) 




^ cosh 



(17) 


dGxx 

~dT~ 


dGp 

dT 


(!-£■)- ao 


dT dOy 


dOy dT 


(23) 
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FIG. 2. Panel (a) The longitudinal conductivity (21) (in 
units of CTo) is plotted as a function of the ratio DvjW for 
different values of the boundary scattering length Ih'- Ih = oo 
(thick solid line), Ih ~ W W (dashed line), Ih = W (dotted 
line), fb = 0.1 IT (dash-dotted line), and Zb = 0 (thin solid 
line). Panel (b) The current-density profile —eJxiy), normal¬ 
ized by the total current I/W, is plotted as a function of y for 
Zb = 0 (no-slip BCs) and different values of D^-. = 0 (thick 

solid line), = 0.01 W (dashed line), and = 0.1 IT (dot¬ 
ted line). The case Di^ ^ IT, corresponding to Poiseuille flow, 
is represented by a thin solid line. 

is the sum of two contributions with opposite signs. The 
first term on the right-hand side of Eq. (23) is negative, 
because < 1 and dao/dT < 0. The latter inequal¬ 
ity holds because the scattering rate describing mo¬ 
mentum non-conserving collisions is a monotonically in¬ 
creasing function of temperature^®. On the contrary, the 
second term is positive, because the dF/dD^, > 0 and 
dD^/dT < 0. The vorticity diffusion length D^, decreases 
with increasing temperature because both v and r are 
decreasing functions of T. We therefore conclude that, 


due to viscosity, axx (Pxx) can increase (decrease) upon 
increasing temperature. This is the so-called Gurzhi ef- 
fect^^. The existence of this effect relies crucially on the 
nature of BCs that are used to solve the hydrodynamic 
equations. In particular, it disappears for free-surface 
BCs. All previous experimental studies of transport in 
graphene and other 2D electron liquids we are aware 
of have reported monotonic temperature dependencies 
(i.e. no evidence of the Gurzhi effect) in the ordinary 
longitudinal geometry in the linear-response regime. We 
therefore conclude that free-surface BCs are the most ap¬ 
propriate for a weak driving current I. In this case, axx 
depends only on the unknown damping rate which 
can therefore be determined from an ordinary four-point 
longitudinal transport measurement at every value of h 
and r, i.e. = e'^nf (maxx)- 

In the next Section, we will discuss another hydro- 
dynamic phenomenon occurring in 2D electron liquids, 
i.e. the formation of whirlpools in electron flow^®, yield¬ 
ing a clear-cut experimental signal of hydrodynamic 
transport in weakly non-local linear-response transport 
measurements. Since the experimental data in Ref. 49 do 
not show any Gurzhi effect in the linear-response regime, 
in the next Section we will utilize only the free-surface 
BCs (Zb —>■ oo). Whirlpools in hydrodynamic electron 
flow, however, do exist also when no-slip BCs are used^®. 
In this sense whirpools are a much more robust phe¬ 
nomenon that the Gurzhi effect in longitudinal trans¬ 
port. Whirpools are also more dramatic in experimental 
appearance. 


IV. NON-LOCAL TRANSPORT AND THE 
IMPACT OF VISCOSITY 

We now present the solution of the problem posed by 
Eqs. (1), (5), (13), and (15) in the rectangle (— 00 , 00 ) x 
[—IT/2, IT/2]. We use free-surface boundary conditions, 
corresponding to Zb —> 00 . 

To this end, we introduce the Fourier transform with 
respect to the longitudinal coordinate x in the three equa¬ 
tions (1) and (5) (the latter for the two components Jx 
and Jy). The Fourier transform of a function f{x, y) will 
be denoted by f{k,y). These equations can be grouped 
into a linear system of three second-order ordinary differ¬ 
ential equations (ODEs) with respect to the independent 
variable y. It is convenient to rewrite this system in terms 
of four first-order ODEs. We find 


/ kJx{k,y) \ 


( ^ 

0 

1 

0 \ 

/ kJx{k,y) \ 

kJy{k,y) 

— h 

—i 

0 

0 

0 ^ 

kJy{k,y) 

dyJx{k,y) 

— r\j 

1 + I/(fcD,)2 

0 

0 

-iUkD^f 

dyJxik, y) 

\ k'^ao(j){k,y)/e / 


V 0 

1 + (kD,)^ 

i{kD,f 

0 / 

\ k'^ao(i>{k,y)/e / 
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Eq. (24) can be solved by diagonalizing the 4x4 ma¬ 
trix on the right-hand side, which has four distinct eigen¬ 
values: ±1 and ±-^1 -I- l/{kD^Y. The general solution 
will therefore be a linear combination of exponentials of 
the form exp(AiA:?/) where Vi (A^) are the eigen¬ 

vectors (eigenvalues) of the matrix. The four unknown 
coefficients can be found by enforcing the desidered 
BCs. These are found by taking the Fourier transform of 
Eqs. (13) and (15) with respect to x: 

Jy{k,y = ±W/ 2 )=J±{k) (25) 

eMk)W 




a—dz 


O-Q 


and 


\OyJx{k,y') 'ikJy{k^yY\y—±^w/2 — 0 


The solution reads as follows: 


w) + w) 


(26) 


(27) 


a=± 


Jl(k,y) = ^c.(k)W X j ik Aa (kw, (kw, 


and 


Jy[k,y) = ^ Mk)W X idy [Aa {kW, + ^Aa {kW, 


cx=± 


VF2 v^'w'w 


+ “^ikPsa, (kW,^, 

1F2 V 'w’wj 


(28) 


(29) 


In writing Eqs. (27)-(29) we have introduced the follow¬ 
ing functions of dimensionless arguments: 


Fi±{k,y) = i 


sinh{ky) ^ cosh{ky) 


fccosh(ft/2) fcsinh(fc/2) 


F2±{k,y) = - 


and 


sinh(fci/) ^ cos\i(ky) 


F^±{k,y,X) = 


cosh(fc/2) sinh(fc/2) 

cosh (y + A“2) 


ik 


cosh(l/2Afc^ + A“2) 


± 


sinh(y + A“ 2 ) 


sinh(l/2Afc^ + A“2) 


(32) 


(30) Eqs. (27)-(32) are the most important results of this Ar¬ 
ticle. 

In general, it is not an easy task to inverse Fourier 
transform Eqs. (27)-(29) to real space, after the func- 

(31) tions J±{k) have been specified. Indeed, this requires to 
calculate a convolution which involves the BCs and the 
functions Fi±{x,y), F 3 ±(x,y), and F' 3 ±(x,y, A) in real 
space. We now introduce the inverse Fourier transforms 
of the functions Fi±{k,y), F 2 ±{k,y), and F’ 3 ±(fc,y, A), 
which read: 


and 


Fi±iS:,y) = ^In 


I- 1 -g -\-2suv{'Ky)e 
I -|- e“27r|a;| _ 2 sin(7ry)e“'^l^l 




1 + e-^’^l^l -h2cos(27rAe"^’^l^l 


2 


F 2 ±(x, y) = I — 7rsin(7ry)e + e |l-t-e — 2[cos(27ry)-|-2]e 

±27re-2^l“l [cos(27rA(l + e"^’"l“l)-h2e-2’^l^ll I x [l-h -h 2cos(27rAe”^’"l“l 


1 



[’‘I 


(33) 


■ (34) 

The functions F 3 ±(x,y, A) do not have simple expressions in terms of elementary functions but can be cast in the 
form of exponentially converging series: 

{ OO _ 

^(2^ + l)(-l)^ cos[(2^ + 
t=o 













































y = t1/2 

\x\ 

y = -l/2 

1 

y = l/2 

\x\ » 1 
y = t1/2 

Fi 

7r-Dn(l=Fe-’^l"l)-f |S|/2 

7r“^ ln(7r|x|) 

7r“^ ln(2) 

\x\/2 =F 

F2- 

±7re-’^l®l/(lTe“"'"l)^ 

{'KX^)~^ 

—ixlA 

±7re-"l*l 

F3- 

(l±l)V(i)/2 

6'{x) 

0 

0 


TABLE I. Explicit expressions of the functions Fm- defined in the main text, evaluated at j/ = =f 1/2- We also summarize 
useful asymptotic behaviors in the limit \x\ <C 1 and |x| S> 1. Similar expressions can also be obtained for the quantities 
Fm+{S;,y) by noting that Fm+{x,y) = Fm-{x,-y). 


OO 

i=i 


(35) 


In the following we will make use of a number of asymp¬ 
totic behaviors of the functions Fi±{x, y), F2±(i, y), and 
F 3 ±{x,y, ^), which have been listed for the sake of con¬ 
venience in Table 1. 

The task of calculating the potential and currents in 
real space simplihes substantially if the currents J7±(x) 
in Eq. (13) can be represented by the sum of a finite 
number of delta-functions in real space. Let us focus on 
the geometry of Fig. 1, where the Fourier transform of 
the BCs (14) reads J_(fc) = - e-*''“o)/e and 

J+{k) = 0. In this case, we find that the steady-state 
current pattern can be written as 


J(r) = —V(/)(r) — nD^V x u:{r) , (36) 

e 


where the potential (f){r) and vorticity ti;(r) = zuj{r) are 
given by: 



and 


lir) = - — 


21 


enW^ 


F {— JL 
' w’ w’ w 


- Fo 


/ x+ 

y 



w’ 

w J 


(38) 


In Eqs. (37) and (38) we have introduce the shorthand 
x± = X ± Xq, with x+ (x_) representing the lateral sep¬ 
aration between the observation point and the collector 
(injector). 

From Eq. (36) we clearly notice an important feature 
of the solution, i.e. for vanishing viscosity the current 
flow is irrotational. More precisely, the viscosity plays 
a twofold role: it modifies the irrotational contribution 
due to the electric potential ())(t) and introduces a finite 


vorticity. It is noteworthy that these effects yield inde¬ 
pendent experimental signatures: the modihcation of the 
electrical potential can be detected by monitoring the re¬ 
sistances in a non-local configuration (or by carrying out 
scanning probe potentiometry), while the vorticity gener¬ 
ates a magnetic field, which can be detected by scanning 
probe magnetometry. These two effects are discussed in 
detail in the following Sections. 


A. Spatial dependence of the 2D electrical 
potential, charge current, and non-local resistances 


Illustrative results for the spatial map of the 2D elec¬ 
trical potential (j){r) —Eq. (37)— and the charge current 
pattern —eJ{r) —Eq. (36) —are shown in Fig. 3. For 
typical values of the drive current I and conductivity, 
i.e. J = 20 /rA through a submicron constriction and 
CTo = 20 mS, we find that the scale over which the 2D 
electrical potential changes is (j)Q = I/uq = 1 mV. We 
clearly see that in the case v ^ 0—panels (b) and (c) in 
Fig. 3— whirlpools with a spatial extension ~ Dy develop 
in the spatial current pattern —eJ(r\ to the right of the 
current injector and to the left of the current collector. 
Once again, the spatial variations of the 2D electrical po¬ 
tential 4>{r) are amenable to experimental studies based 
on scanning probe potentiometry. 

In passing, we note that near the current injector at 
X = Xq the potential is dominated by the singular parts 
of the functions Fm-- Taking the limit W/xq —>■ oo in 
Eq. (37) we find an extremely simple expression for the 
potential near the injector: 




I 

TTCTo 





(39) 


where r is the distance from the injection point, 6 is 
the angle measured from the injection direction y, and 
i? is a length determined by BCs far from the contact. 
Note that changing R is equivalent to changing (/) by an 
arbitrary additive constant. 

If dc transport is to be used as the main tool to detect 
hydrodynamic electron flow, it is pivotal to understand 





















the spatial dependence of the non-local resistance i?NL, 
which we define in the following way: 
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p / N _ y) - (/'(a: + 00 , y) 

Rni.{x,y) = - z - 


(40) 


where the quantity (j)(x —> +oo,y) does not depend on 
y. Because of (37), we find that, at each point in space, 
R^ij[x,y) is a quadratic function of D^-. 


R^'Li.x, y)ao = a[x, y)Dl + b{x, y) , (41) 


where 






and 




To make contact with Ref. 49, we now introduce the 
“vicinity” resistance, which is the non-local resistance 
measured on the edge where current is injected, at a dis¬ 
tance Ax from the current injector: 


i?v(Ax) = RNL(a;o + Aa;, -W/2) . (44) 


Using Eqs. (41), (42), (43), the asymptotic results in Ta¬ 
ble 1, and taking the limit Xq lU, we find 


i?v(Ax) = 


-27re“’"A^I/w^ 
] 4/2 ('x _ g-7i-|Aa;|/W) 


-I- 

2 ^ 


+ 


In (l - + Ax 0(-Ax) 


(45) 


Here 0(x) is the Heaviside step function. For positive Ax 
the two terms in Eq. (45) have opposite sign. For this 
reason, Ry{Ax) is expected to change sign as a function 
of D,y. The change of sign of the vicinity resistance is a 
key signature of the viscous contribution to the electric 
potential. Maximum sensitivity to viscosity is achieved 
when slightly non-local or “vicinity” voltage drops are 
measured outside the region [—xo;xo] where the current 
flux is maximum. The vicinity resistance (45) rapidly 
decays for |Ax| 3> TU/tt: it is therefore pivotal to mea¬ 
sure^® the potential (j>{x, y) for a lateral separation Ax 
from the current injection point which is of the order of 
the vorticity diffusion length D^. 

Once the rate of momentum-non-conserving collisions 
T~^(n,T) is measured from an ordinary four-point lon¬ 
gitudinal measurement of Gxx (as explained in Sect. HI), 
a measurement of the vicinity resistance Ry{Ax) yields 
a map v = v{n, T) of the kinematic viscosity of the 
2D electron liquid. We hasten to stress that our all¬ 
electrical non-local protocol to measure the kinematic 
viscosity v = v(n, T) applies to any 2D electron liquid 
driven into the hydrodynamic transport regime (and not 
only to doped SLG and BLG). 



FIG. 3. (Color online) Steady-state spatial map of the 2D 
electrical potential (j>{r) (in units of 4>o = 7/co) and charge 
current streamlines —eJ{r) in a Hall bar device like the one 
depicted in Fig. 1 with xo = W. Different panels refer to 
different values of the vorticity diffusion length Di, = 0 
[panel (a)], = 0.5 W [panel (b)[, and = IT [panel 

(c)]. Whirlpools are clearly seen in the bottom right and 
bottom left of panels (b) and (c). No whirlpools occur in the 
absence of viscosity, as in panel (a). In each panel, the current 
streamlines change color from white (high current density) to 
black (low current density). 


For the sake of completeness, we note that the authors 
of Ref. 24 have proposed an ac Gorbino disk viscometer, 
which allows a determination of the hydrodynamic shear 
viscosity from the dc potential difference that arises be¬ 
tween the inner and the outer edge of the disk in response 
to an oscillating magnetic flux. 


B. Spatial depedence of the current-induced 
magnetic field 

Because the steady-state current —eJ{r) generates a 
magnetic field in the proximity of 2D electron system, 
whirlpools and viscosity-dominated hydrodynamic trans¬ 
port can also be detected by scanning probe magnetom- 
etry (see, for example. Refs. 58-60). 

As shown in Appendix A, in a sample in which a back- 
gate is placed at a distance z = —d below the graphene 
sheet with d ^ W, D^, a local relation exists between the 
z-component Bz{r,z > 0) of the magnetic field and the 
vorticity uj{r). In SI units, this relation reads as following 

Bz{r,z > 0) = —e/iohdaj(r) , (46) 

where /tq = 47r x 10“^ N/A® is the magnetic constant 
and the vorticity has been introduced earlier in Eq. (38). 

A typical 2D spatial map of Bz{r,z > 0) is shown in 
Fig. 4 for different values of the vorticity diffusion con¬ 
stant Djy. In this figure the magnetic field is plotted in 
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FIG. 4. (Color online) Spatial map of the £-component of 
the magnetic field, Bz{r,z) (in units of Bo = fioId/W^), 
generated by the 2D steady-state current pattern J{r) and 
calculated immediately above the graphene sheet, i.e. for 
0 < 2 <C W,Dv These results have been obtained for the 
same parameters as in Fig. 3. Different panels refer to differ¬ 
ent values of D^-. = 0.5 W [panel (a)] and = W 

[panel (b)[. In this figure we have not shown results for 
= 0: in the absence of viscosity B^{r,z) is identically 
zero. 

units of Bq = ^old/W^. For / = 200 fiA, W = 1 /im, 
and d = 80 nm, we find Bq = 20 /iT. This value is well 
within reach of current technology^®^®*^. 

V. SUMMARY AND FUTURE PERSPECTIVES 

In summary, we have presented a theoretical study of 
dc transport in graphene driven into the hydrodynamic 
regime. As highlighted in Ref. 49, this regime occurs in 
a wide range of temperatures and carrier concentrations. 

Our theory, which applies only to the doped regime, re¬ 
lies on the continuity (1) and Navier-Stokes (5) equations, 
augmented by suitable boundary conditions at Hall bar 
edges. We have demonstrated analytically that a combi¬ 
nation of ordinary four-point longitudinal transport mea¬ 
surements and measurements of non-local resistances in 
Hall bar devices can be used to extract the hydrodynamic 
shear viscosity of the two-dimensional electron liquid in 
graphene"^®. 

We have also discussed how to probe the viscosity- 
dominated hydrodynamic transport regime by scanning 
probe methods. Indeed, we believe that it possible to 
observe hydrodynamic electron flow with spatial reso¬ 
lution by using available scanning probe potentiometry 
and magnetometry setups. Spatial maps of the two- 
dimensional electrical potential ^’{r) and current-induced 
magnetic field Bz(r,z > 0) for experimentally relevant 
parameters are shown in Figs. 3 and 4. 

We wish to emphasize that our theoretical approach is 
immediately applicable to any 2D electron liquid in the 
hydrodynamic transport regime. 

In the future, we plan to extend our theoretical ap¬ 
proach to semimetals with linear and quadratic band 
touchings at the charge neutrality point, where viscos¬ 


ity is expected to be very low^^ and the Reynolds num¬ 
ber (16) is expected to be very large, possibly enabling 
the observation of electronic turbulence. This will require 
to deal with thermally-excited carriers, coupling between 
charge and heat flow, and non-linear terms in the Navier- 
Stokes equation. We also would like to gain a fully micro¬ 
scopic understanding of momentum-non-conserving col¬ 
lisions in the hydrodynamic transport regime by treat¬ 
ing smooth scalar and vector potentials due to disorder 
(strain®^, charged impurities®®, etc.) along the lines of 
what was done by the authors of Ref. 19 for the case of 
a smooth scalar potential. 

Last but not least, we strongly believe that hydrody¬ 
namic flow and the shear viscosity of 2D electron liquids 
can also be accessed®^ by scattering-type near-held opti¬ 
cal spectroscopy (see, for example. Ref. 38 and references 
therein to earlier work) in the Terahertz spectral range, 
since this technique measures the non-local conductiv¬ 
ity a{q,u}). Terahertz radiation is required a) to make 
sure that the hydrodynamic inequality ujTee <C 1 is satis- 
hed (with Tee = ^ee/vp) and b) to have measurable non¬ 
local effects due to viscosity, since the latter decreases 
quickly^® as a function of the external probe frequency 

UJ. 


Appendix A: Local inductance approximation 

In this Appendix we present a derivation of Eq. (46). 
We use SI units and the Coulomb gauge V • A = 0 for 
the vector potential. We assume that a bottom gate po¬ 
sitioned at 2 = —d exists below the graphene sheet. This 
will be treated as a perfect conductor. 

The vector potential is related to the steady-state cur¬ 
rent pattern by the following 3D Poisson equation: 

(v^ + A(r, z) = /roeJ ir)S{z) . (Al) 

This is similar to the Poisson equation in Eq. (12) for the 
scalar potential <I>(r, 2 ). 

Fourier transforming Eq. (Al) with respect to the in¬ 
plane coordinate r we find: 

(- 9 ^ + ^)-4(q,2) =/ioeJ(q)d( 2 ) , (A2) 

where q ■ J{q) = 0 because of the continuity equation 
Eq. (1). 

The general solution of Eq. (A2) is: 

-4(9,2) = -/ioeJ(9)^^^^ +d+(q)e«^ -f a_(q)e"«^ . 

(A3) 

The quantities a±{q) must obey the condition ^± 2 ( 9 ) = 
T *9 • ^±{q)/q to enforce the Coulomb gauge and must 
be determined from BCs. Requiring A( 9 , 2 —)■ -boo) = 0 
implies that d+( 9 ) must vanish. 
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The corresponding z component of the magnetic field 
is given by: 

B^{q,z) = [iq X A{q,z)]^ . (A4) 

Since must vanish at the gate position, i.e. for z = —d, 
we find: a^{q) = fioej{q)e~‘^'>‘^/{2q). In deriving the 
previous result we have assumed that a^{q) ■ q = Q. 

The Fourier transform of the vector potential is there¬ 
fore given by: 

^ , p-9bl _ p-qi'^d+z) 

= -noeJ{q) - — - . (A5) 

Now, if d and \z\ are small with respect to the lateral 
lengthscales W and Dp over which the steady-state cur¬ 
rent pattern J{r) changes in the sample, i.e. d,\z\ ^ 
W, Dp , the above formula can be approximated for z > 0 
as: 

A{q,z > 0) « —e/xodJ(q) . (A6) 


Making use of Eq. (A4) and transforming back to real 
space we finally obtain the desired result. 


Bz{r, z > 0) « —e/xod[V xJ{r)]z = —e^Qnduj{r) . (A7) 
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